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In this paper a method of analysis of dynamic response of laminated 
composite plates is presented. The analysis is carried by using a 
hybrid- stress finite -element numerical technique established by the 
authors in their earlier publication. By using this approach the 
response of simply supported laminated plates subject to sinusoidal 
loading are investigated. For the solution of the finite -element 
equations of motion of free vibrations and dynamic response problems , 
two effective methods of solution, the space iteration method and the 
Newmark direct integration method are used. These two methods are 
discussed in this paper. 

INTRODUCTION 

Since Pian [1] first established the assumed stress hybrid finite element 
model and derived the corresponding element stiffness matrix in 1964, the hybrid 
stress model has been shown highly accurate, and easy to fulfill the 
compatibility condition of the finite element method. Laminated thick plate 
element has been developed by Mau et al. [2] by using hybrid stress method. In 
the comparison of results with elasticity solution [3,4], they observed excellent 
accuracy in predicting both displacements and stresses. In their assumption for 
the stress field, transverse normal stress was not included. Constant transverse 
displacement through the laminate thickness was also assumed. These assumptions 
did not agree well with the actual mechanism of deformation of laminated plates 
in bending. Spilker [5] developed an eight-node isoparametric multilayer plate 
element for the analysis of thin to thick fiber- reinforced composite plates. 
This model has the generality in describing laminate response and can be easily 
used to implement to attack complex laminated plate problems, but the assumption 
of constant transverse displacement through laminate thickness still remains. 

The hybrid stress model is based on the modified complementary energy 
principle. An optimum choice of the number of the assumed stress modes for given 
boundary displacement approximation can be made, which give greater flexibility 
in the descriptions of the stress field. The detail of the development of this 
method is documented in [6]. 

In the present investigation, a three-dimensional eight-node hybrid stress 
element has been developed to analyze free vibrations of laminated plates. All 
six stress components are included and assumed independently within each layer 
through stress polynomials with 55 unknown stress parameters. The stress field 
within each layer satisfies the dynamic equilibrium equations of free vibration. 
The interface traction continuity and laminate upper/lower surface traction- free 
conditions are also enforced. The displacement field is interpolated in terms of 
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nodal displacements through shape functions. The displacements (u,v,w) are 
assumed to vary linearly through the thickness of each lamina. 

To solve the governing finite element equations of motion for a linear 
dynamic analysis without damping the well-known Newmark direct integration method 
[7] will be used to integrate the following equation 

[M]{V) + [k] { q } - [Q] 

step by step. 

The dynamic response of simply supported laminated composite plates under a 
dynamic sinusoidally distributed load 

Q - Q sin ^ sin ~^H(t) 
o a b 

is analyzed. Numerical results of [90/0] antisymmetric cross-ply laminate and 
[0/90/0] symmetric cross-ply laminate are presented. Center deflection, bending 
stresses a x , ay, transverse shear stresses r xz and r y Z and normal stress a z for 
both laminates are plotted as a function of time. Fast convergence is observed . 

SUBSPACE ITERATION METHOD 

the free vibration finite element equations of motion with damping neglected 

are 


M q* + K q - 0 (1) 

where K is the structure stiffness matrix and M is the structure mass matrix. 
Equation {1} can be solved by expressing the field variable as 

q - *e iwt (2) 

where <f> is a nodal vector of order n, t the time variable, and u> the natural 

frequency of vibration of the plate in the mode described by the vector 

<f > . Substituting equation (2) into Equation (1) yields the generalized eigenvalue 

problem 


K <f> - w 2 M <j> - 0 (3) 

from which </> and c*> can be determined. For matrices of dimension n x n, there 

will be n eigensolutions , <f>\) , (o>2^ , <f> 2 ) , (w n ^ , ^ n ) * An important 

property of the eigenvectors is that they satisfy the orthogonality conditions, 
i.e. 


<t> i T M - S ij 

(4) 

K $ij 

and 0 < < 

There are many different techniques existing for the solution of eigenvalue 
problems. Since the procedures for the eigenvalues problems are time consuming, 
the choice of an appropriate and effective method is an important factor for the 
general application, especially in the large eigenvalue problem. The subspace 
iteration method suggested by Bathe [8] will be adopted to conduct the 
investigation. This method has been used extensively in a number of general- 
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purpose finite element analysis programs and has proven cost-effective and 
reliable. In structure analysis, the lowest few eigenvalues (natural 
frequencies) are the main concern of investigators. The basic objective in the 
subspace iteration method is to solve for the p smallest eigenvalues and 
corresponding eigenvectors, which satisfy 

K <f> — M <f> A (5^ 

where <f> - <f> 2 , <£p] 

and A is a diagonal matrix of an< j the eigenvector 4>t also satisfies the 

orthogonality conditions (Equation (4)). 

The subspace iteration method consists of three steps [8]: 

1. Establish q starting iteration vectors; q > p, q - min(2p, p+8) is a 
ptoper selection, where p is the number of eigenvalues and eigenvectors 
to be calculated. 

2. Use simultaneous inverse iteration on the q vectors and Ritz analysis 
to extract the "best" eigenvalue and eigenvector approximations from 
the q iteration vectors . 


For k - 1, 2 

K x k+l - M X k (6) 

where is the starting iteration vector. 

Find the projections of the operators K and M, 

K k+1 “ W K X k+1 

(7) 

M k+1 “ W M X k+1 

Solve for the eigensystem of the projected operators, 


K k+1 Q k +1 “ M k+1 Q k+1 A k+1 (g) 

Find an improved approximation to the eigenvectors, 

K k+1 “ x k+l Qk+1 (9) 

As k -*■ <f> 

A k+1 - A and X k+1 -+ <j> (10) 

3. After iteration convergence, use the Sturm sequence check to verify 
that the required eigenvalues and corresponding eigenvectors have been 
calculated. 

Reference [8] has presented very detailed descriptions about the subspace 

iteration method and is a good reference to use to become familiar with this 
method. 
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NEWMARK DIRECT INTEGRATION METHOD 


The governing finite element equations of motion for a linear dynamic 
analysis are 


Mq + Kq»0 < n > 

where M and K are the mass and stiffness matrices, q and q are the nodal 
displacement and acceleration vector, and Q is the load (nodal force) vector of 
the finite element system. In finite element analysis, there exist many 

effective numerical procedures to solve the linear differential equations, 
Equation (11). Basically, they can be divided into two methods of solution: the 
direct integration method and the mode superposition method. In the present 
study, the Newmark direct integration method [7] will be followed to integrate 
Equation (11) step by step. The following assumptions are made in the numerical 
analysis : 

q n +l “ 5n + [ (l-5)Vn +5 qn+l] At 

q n +i - q n + q n At + K \ -«)Vn + “qn+i] At2 (13) 

where At is the time step size, n is the step number, and the parameters S and a 
control integration accuracy and stability. At time t n+ i - (n+l)At, the finite 
element equations of motion (Equation (11)) are described as: 


M qn+1 + K q n +i ” Qn+1 


(14) 


Solving from Equation (13) for q^+l in terms of q n+1 , and then substituting into 
Equation (14) and rearranging the terms transforms the equations to the form 


A A 

K qn+1 ” Qn+1 


(15) 


where 


A 

K - a Q M + K 


(16) 


Qn+1 ” Qn+1 + M(a 0 q n +aiq n +a2q n ) 


and 

a Q - l/a(At) 2 
a.\ - 1/a (At) 


(17) 


a 2 — (1 - 2a)/2a 


Once the displacements q n +i at time step n+1 are known from Equation (15) , the 
velocities and accelerations can be computed by using Equations, (12) and (13) 
and expresses as 


• t 

q n +i - 


a o ( qn+i 


q n ) - a l5n - a 2 *qn 


(18) 


q n +i - q n + a 3qn + a 4q’ n +i 


(19) 


where 
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a 3 - (1 - S) (At) 
a4 - 6 (At) 


( 20 ) 


A special scheme originated by Newmark with S = 0.5 and a - 0.25 is used here to 
integrate the equations step by step. These values correspond to the constant- 
average -acceleration method, which gives an unconditionally stable numerical 
scheme [7]. 

DYNAMIC RESPONSE OF A SIMPLY SUPPORTED LAMINATED SQUARE PLATE 

The dynamic response of simply supported laminated plates is presented in 
this section. The laminates are subjected to suddenly applied sinusoidally 
distributed pulse loading, 

q(x, y, t) = (q 0 sin7rx/a sin?ry/b)H( t) (21) 

where H(t) is the Heavyside step function. The following two laminated plates 
are considered: 

1. A two-layer anti -symmetric cross-ply (90/0) square laminate with layers 
of equal thickness. 

2. A three-layer symmetric cross-ply (0/90/0) square laminate with layers 
of equal thickness. 

In both problems, the same material properties as in Putcha and Reddy [9] are 
employed for each individual layer. 

E l = 525 GPa 

E t = 21 GPa 

G lt - G tt =10.5 GPa (22) 


V LT “ V TT “ 0.25 
p = 0.8 g/cm^ 

Owing to the biaxial symmetry of the laminate geometry, only one quadrant of the 
laminate is analyzed. The geometry configurations and boundary conditions of the 
finite element model are shown in Figure 1. 
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(a) Top view and boundary conditions of laminate 
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(b) Two-layer anti-symmetric laminate 
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(c) Three-layer symmetric laminate 


Figure 1 Laminate geometry configurations and boundary 
conditions 


The normalized deflection and stresses are described as 
w = 1000E-ph^w/q o a^ 

(ct x , CTy) = 10(<7 X , a y )/q Q S 2 (23) 

( a z » r xz , r yz) = r xz » r yz)/ c lo^ 

S * a/h z - z/h 

In the present study q Q is taken to be 100 and the time step size is equal to 5 
microseconds. The normalized central transverse deflection, w(a/2, a/2, 0), as a 
function of time for a two-layer anti-symmetric simply supported cross-ply square 
laminate under sinusoidal loading is shown in Figure 2. Throughout Figures 3, 4, 
5 and 6, the stresses with respect to time for a two- layer laminate are plotted. 
In Figures 3 and 4, it is observed that the normalized normal stress a x at center 
top surface of the laminate is close to the normalized normal stress Oy at center 
bottom surface of the laminate, except a x is in tension and <jy is in compression. 
As shown in Figure 5, the variation of normalized shear stress 7 XZ is similar to 
7y Z . From the plots, it is seen that the periods of the transient response for 
w, a x , ay, r xz , and 7y Z are very closely related. This fact agrees with the 
results of Putcha and Reddy [9] . The period for the normalized transverse normal 
stress a z (a/2, a/2, h/2) , is much shorter when compared with others. A shorter 
time step size (At = 1 microsecond) is employed to observe the periodic response 
of the normalized transverse normal stress a z as shown in Figure 6. 
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Figure 5 


Transverse shear stresses versus time for a 2-layer (90/0) laminate 
under suddenly applied sinusoidal loading 



Figure 6 Transverse normal stress versus time for a 2-layer (90/0) laminate 
under suddenly applied sinusoidal loading 

The normalized central transverse deflection with respect to time for a 
three -layer simply supported cross -ply (0/90/0) square laminate under sinusoidal 
loading is shown in Figure 7. Figure 8 contains the normalized maximum normal 
stresses, ct x and Oy, as a function of time. The normalized shear stresses, r xz 
and r yz , are shown in Figure 9. The periods for the normalized deflection and 
stresses are similar. In Figure 9, the amplitude of the response is larger for 
r xz than for r yz ; it is because the bending stiffness is higher in the x- 
direction than in the y- direction for a three -layer (0/90/0) laminate with layers 
of equal thickness. The normalized central normal stress distribution, a x , 

through the thickness of the laminate for time from 20 to 80 microseconds is 
shown in Figure 10. 

In the present study, a four-node isoparametric plate element with 48 
assumed stress parameters for each lamina is used. Fast convergence is observed; 
only a 5 x 5 mesh is modeled in a quadrant of the laminate. 
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Center deflection versus time for a 3-layer cross-ply (0/90/0) 
laminate under suddenly applied sinusoidal loading 
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Normal stresses versus time for a 3-layer (0/90/0) laminate under 
suddenly applied sinusoidal loading 
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Figure 9 Transverse shear stresses versus time for a 3-layer (0/90/0) laminate 

under suddenly applied sinusoidal loading 
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Figure 10 Through- thickness center normal stress versus z for a 3-layer 
(0/90/0) laminate under suddenly applied sinusoidal loading 
for tine from 20 to 80 microseconds 


CONCLUDING REMARKS 


The proposed three-dimensional hybrid stress finite element method in 
conjunction with the Newmark's direct integration method seems to be a powerful 
technique for analyzing laminated composite plates under dynamic loading. By 
using this approach the transverse deflection, the in-plane bending stresses and 
the interlaminar shear stresses and normal stress can be evaluated very easily 
without consuming too much computation time. This method can also be used to 
analyze forced vibration problems of laminated composite under impact loading. 
The result will be published in the near future. 


484 


REFERENCES 


1. Pian, T. H. H. , "Derivation of Element Stiffness Matrices by Assumed Stress 
Distribution," AIAA J . Vol . 2, pp. 1333-1336, 1964. 

2. Mau, S. T., Tong, P., and Pian, T. H. H. , "Finite Element Solutions for 
Laminated Thick Plates," J. Composite Materials . Vol. 6, pp . 304-311. 

3. Reissner, E. , and Stavsky, Y. , "Bending and Stretching of Certain Types of 
Heterogeneous Allotropic Elastic Plates," J. Applied Mechanics . Vol. 28, pp. 
402-408, 1961. 

4. Whitney, J. M. , and Panago , N. J., "Shear Deformation in Heterogeneous 
Anisotropic Plates," J. Applied Mechanics Transactions, ASME, 92, pp. 1031- 
1036, 1970. 

5. Spilker, R. L. , , "A Hybrid-Stress Finite Element Formulation for Thick 
Multilayer Laminates," Computer and Structures . Vol. 11, pp . 507-514, 1980. 

6. Sun, C. T. , and Liou, W. J., "A Three-Dimensional Hybrid Stress 
Isoparametric Element for the Analysis of Laminated Composite Plates," 
Computer and Structures Vol. 25, No. 2, pp. 241-249, 1987. 

7. Newmark, N. M. , "A Method of Computation for Structural Dynamics," ASCE 
J. Engineering Mechanics Division . 85, pp. 67-94, 1959. 

8. Balke, K. J., Finite Element Procedures in Engineering Analysis . Prentice - 
Hall, Englewood Cliffs, New Jersey, 1982. 

9. Putcha, N. S., and Reddy, J. N. , "On Dynamics of Laminated Anisotropic 
Plates Using a Refined Mixed Plate Element" ASME Winter Annual Meeting, New 
Orleans, La. pp. 161-169, 1984. 


485 


